rm(list = ls())
load("dataBJPOLS.RData")

inst.1 <- "judgeiv_hd * as.factor(race4) * as.factor(judge_exp_cat)"
endo.1 <- "pti * as.factor(race4) * as.factor(judge_exp_cat)"
outc.1 <- "vote_post"

time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(severity)"
demo.controls <- "age + I(age^2) +  as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before"
case.controls <- "as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)"

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

resultsB <- resultsH <- resultsW <- resultsM <- resultsHM <- resultsWM <- list()
boot.run <- T

if(boot.run) {
   set.seed(1231)
   for(i in 1:500) {
      out <- last.cases[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
      last.cases.boot <- last.cases[last.cases$this_id_fa %in% out$V1, ]
      
      m1a1 <- ivreg(form.1, data = last.cases.boot)
      m1a2 <- ivreg(form.2, data = last.cases.boot)
      m1a3 <- ivreg(form.3, data = last.cases.boot)
      
      resultsB[[i]] <- c(m1a1$coefficients["pti"], 
                         m1a2$coefficients["pti"], 
                         m1a3$coefficients["pti"])
      
      resultsW[[i]] <- c(m1a1$coefficients["pti:as.factor(race4)W"], 
                         m1a2$coefficients["pti:as.factor(race4)W"], 
                         m1a3$coefficients["pti:as.factor(race4)W"])
      
      resultsH[[i]] <- c(m1a1$coefficients["pti:as.factor(race4)H"], 
                         m1a2$coefficients["pti:as.factor(race4)H"], 
                         m1a3$coefficients["pti:as.factor(race4)H"])
      
      resultsM[[i]] <- c(m1a1$coefficients["pti:as.factor(judge_exp_cat)1"], 
                         m1a2$coefficients["pti:as.factor(judge_exp_cat)1"], 
                         m1a3$coefficients["pti:as.factor(judge_exp_cat)1"])
      
      resultsHM[[i]] <- c(m1a1$coefficients["pti:as.factor(race4)H:as.factor(judge_exp_cat)1"], 
                          m1a2$coefficients["pti:as.factor(race4)H:as.factor(judge_exp_cat)1"], 
                          m1a3$coefficients["pti:as.factor(race4)H:as.factor(judge_exp_cat)1"])
      
      resultsWM[[i]] <- c(m1a1$coefficients["pti:as.factor(race4)W:as.factor(judge_exp_cat)1"], 
                          m1a2$coefficients["pti:as.factor(race4)W:as.factor(judge_exp_cat)1"], 
                          m1a3$coefficients["pti:as.factor(race4)W:as.factor(judge_exp_cat)1"])
   }
   
   SEsB <- round(apply(do.call('rbind', resultsB), 2, sd), 3)
   SEsW <- round(apply(do.call('rbind', resultsW), 2, sd), 3)
   SEsH <- round(apply(do.call('rbind', resultsH), 2, sd), 3)
   SEsM <- round(apply(do.call('rbind', resultsM), 2, sd), 3)
   SEsWM <- round(apply(do.call('rbind', resultsWM), 2, sd), 3)
   SEsHM <- round(apply(do.call('rbind', resultsHM), 2, sd), 3)
   
   SEsBMa <- round(apply(do.call('rbind', resultsB) + do.call('rbind', resultsM), 2, sd), 3)
   SEsBNo <- round(apply(do.call('rbind', resultsB), 2, sd), 3)

   SEsWMa <- round(apply(do.call('rbind', resultsB) + do.call('rbind', resultsM) + do.call('rbind', resultsW) + do.call('rbind', resultsWM), 2, sd), 3)
   SEsWNo <- round(apply(do.call('rbind', resultsB) + do.call('rbind', resultsW), 2, sd), 3)
   
   SEsHMa <- round(apply(do.call('rbind', resultsB) + do.call('rbind', resultsM) + do.call('rbind', resultsH) + do.call('rbind', resultsHM), 2, sd), 3)
   SEsHNo <- round(apply(do.call('rbind', resultsB) + do.call('rbind', resultsH), 2, sd), 3)
   
   SEsWBMa <- round(apply(-1 * (do.call('rbind', resultsB) + do.call('rbind', resultsM) + do.call('rbind', resultsW) + do.call('rbind', resultsWM))
                          + (do.call('rbind', resultsB) + do.call('rbind', resultsM)), 2, sd), 3)
   SEsWBNo <- round(apply(-1 * (do.call('rbind', resultsB) + do.call('rbind', resultsW))
                          + (do.call('rbind', resultsB)), 2, sd), 3)

   SEsHBMa <- round(apply(-1 * (do.call('rbind', resultsB) + do.call('rbind', resultsM) + do.call('rbind', resultsH) + do.call('rbind', resultsHM))
                          + (do.call('rbind', resultsB) + do.call('rbind', resultsM)), 2, sd), 3)
   SEsHBNo <- round(apply(-1 * (do.call('rbind', resultsB) + do.call('rbind', resultsH))
                          + (do.call('rbind', resultsB)), 2, sd), 3)
   
   SEsHWMa <- round(apply(-1 * (do.call('rbind', resultsB) + do.call('rbind', resultsM) + do.call('rbind', resultsH) + do.call('rbind', resultsHM))
                          + (do.call('rbind', resultsB) + do.call('rbind', resultsM) + do.call('rbind', resultsW) + do.call('rbind', resultsWM)), 2, sd), 3)
   SEsHWNo <- round(apply(-1 * (do.call('rbind', resultsB) + do.call('rbind', resultsH))
                          + (do.call('rbind', resultsB) + do.call('rbind', resultsW)), 2, sd), 3)
}

m1a1 <- ivreg(form.1, data = last.cases)
m1a2 <- ivreg(form.2, data = last.cases)
m1a3 <- ivreg(form.3, data = last.cases)

m1a1_d <- summary(m1a1, vcov = sandwich, diagnostic = T)
m1a2_d <- summary(m1a2, vcov = sandwich, diagnostic = T)
m1a3_d <- summary(m1a3, vcov = sandwich, diagnostic = T)

res1_pti <- coefficients(m1a1)[grep("pti", names(coefficients(m1a1)))]
res2_pti <- coefficients(m1a2)[grep("pti", names(coefficients(m1a2)))]
res3_pti <- coefficients(m1a3)[grep("pti", names(coefficients(m1a3)))]

## White-Black Defendant
c1_wo <- c(-1 * (res3_pti["pti:as.factor(race4)W"]), SEsWBNo[3])
c1_wm <- c(-1 * (res3_pti["pti:as.factor(race4)W"] + res3_pti["pti:as.factor(race4)W:as.factor(judge_exp_cat)1"]), SEsWBMa[3])

## Hispanic-Black Defendant
c1_ho <- c(-1 * (res3_pti["pti:as.factor(race4)H"]), SEsHBNo[3])
c1_hm <- c(-1 * (res3_pti["pti:as.factor(race4)H"] + res3_pti["pti:as.factor(race4)H:as.factor(judge_exp_cat)1"]), SEsHBMa[3])

## Hispanic-White Defendant
c1_hwo <- c(-1 * (res3_pti["pti:as.factor(race4)W"] - res3_pti["pti:as.factor(race4)H"]), SEsHWNo[3])
c1_hwm <- c(-1 * (res3_pti["pti:as.factor(race4)W"] + res3_pti["pti:as.factor(race4)W:as.factor(judge_exp_cat)1"] - 
              (res3_pti["pti:as.factor(race4)H"] + res3_pti["pti:as.factor(race4)H:as.factor(judge_exp_cat)1"])), SEsHWMa[3])

m1 <- data.frame(rbind(c1_wm,
                       c1_wo))

size <- 20
colnames(m1) <- c("estimate", "std.error")
m1$term <- c("t1", "t2")
m1_df <- data.frame(as.matrix(m1)) 
m1_df$estimate <- as.numeric(as.character(m1_df$estimate))
m1_df$std.error <- as.numeric(as.character(m1_df$std.error))

g3 <- {dwplot(m1_df, 
              dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
              vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
      relabel_predictors(c(t1 = "Experienced", t2 = "Inexperienced")) +
      theme_classic() + xlab("Difference in the Effect of\n Pretrial Incarceration on\n Turnout (Black - White)") + ylab("Judge's Experience") + xlim(-.4, .2)+
      ggtitle("Black vs White\nDefendants") +
      theme(plot.title = element_text(face="bold"), legend.position = "none",
            axis.text.x  = element_text(size = size),
            axis.text.y  = element_text(size = size),
            text = element_text(size = size))} 

m1 <- data.frame(rbind(c1_hm,
                       c1_ho))

colnames(m1) <- c("estimate", "std.error")
m1$term <- c("t1", "t2")
m1_df <- data.frame(as.matrix(m1)) 
m1_df$estimate <- as.numeric(as.character(m1_df$estimate))
m1_df$std.error <- as.numeric(as.character(m1_df$std.error))

g1 <- {dwplot(m1_df, 
              dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
              vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
      relabel_predictors(c(t1 = "Experienced", t2 = "Inexperienced")) +
      theme_classic() + xlab("Difference in the Effect of\n Pretrial Incarceration on\n Turnout (Black - Hispanic)") + ylab("Judge's Experience") + xlim(-.4, .2)+
      ggtitle("Black vs Hispanic\nDefendants") +
      theme(plot.title = element_text(face="bold"), legend.position = "none",
            axis.text.x  = element_text(size = size),
            axis.text.y  = element_text(size = size),
            text = element_text(size = size))} 


m1 <- data.frame(rbind(c1_hwm,
                       c1_hwo))

colnames(m1) <- c("estimate", "std.error")
m1$term <- c("t1", "t2")
m1_df <- data.frame(as.matrix(m1)) 
m1_df$estimate <- as.numeric(as.character(m1_df$estimate))
m1_df$std.error <- as.numeric(as.character(m1_df$std.error))

g4 <- {dwplot(m1_df, 
              dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
              vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
      relabel_predictors(c(t1 = "Experienced", t2 = "Inexperienced")) +
      theme_classic() + xlab("Difference in the Effect of\n Pretrial Incarceration on\n Turnout (Hispanic - White)") + ylab("Judge's Experience") + xlim(-.4, .2)+
      ggtitle("Hispanic vs White\nDefendants") +
      theme(plot.title = element_text(face="bold"), legend.position = "none",
            axis.text.x  = element_text(size = size),
            axis.text.y  = element_text(size = size),
            text = element_text(size = size))} 

pdf(file = "./Figures/Figure_05_BJPOLS.pdf", width = 18, height = 6)
ggarrange(g3, g4, g1, ncol = 3, nrow = 1)
dev.off()

